assessing impact of analysis approaches: all_retest
## libraries
library(colorout)
library(here)
library(data.table)
source(here("src", "stroop-rsa-pc.R"))
#render_report("optim_analysis", src_dir = "group_analyses/optim_2021-11-28", params = list(subjlist = "ispc_retest"))
#render_report("optim_analysis", src_dir = "group_analyses/optim_2021-11-28", params = list(subjlist = "all_retest"))
library(knitr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggridges)
library(abind)
library(httpgd)
## settings
hgd()
opts_chunk$set(echo = TRUE)
theme_set(theme_bw(base_size = 10))
## constants
#subjlist <- "all_retest"
subjlist <- params$subjlist
if (subjlist == "ispc_retest") {
session <- "reactive"
} else if (subjlist == "all_retest") {
session <- "proactive"
}
prewhs <- c("none", "obsall", "obsresampbias")
glms <- c("lsall_1rpm", "lssep_1rpm", "condition_1rpm")
roiset <- "Schaefer2018Dev"
measures <- c("cveuc", "crcor")
subjects <- fread(here("out", paste0("subjlist_", subjlist, ".txt")))[[1]]
atlas <- read_atlas(roiset)
## read subject RDMs (just lsall_1rpm):
rdm <- enlist(combo_paste(measures, prewhs, sep = "__"))
for (prewh in prewhs) {
for (measure in measures) {
rdm[[paste0(measure, "__", prewh)]] <- read_rdms(
.measure = measure, .prewh = prewh,
.glmname = "lsall_1rpm", .roiset = roiset, .waves = waves,
.session = session, .subjects = subjects,
.ttypes1 = intersect(ttypes_by_run[[session]]$run1, ttypes[["bias"]]),
.ttypes2 = intersect(ttypes_by_run[[session]]$run2, ttypes[["bias"]]),
.rois = names(atlas$rois)
)
}
}
## read regression weights:
dat <- enlist(combo_paste(measures, glms, sep = "__"))
for (measure in measures) {
for (glmname in glms) {
for (prewh in prewhs) {
if (prewh != "none" && glmname != "lsall_1rpm") next
if (measure == "crcor" && glmname == "condition_1rpm") next
file_name <-
construct_filename_weights(
measure = measure, subjlist = subjlist, glmname = glmname, roiset = roiset, prewh = prewh
)
dat[[paste0(measure, "__", glmname, "__", prewh)]] <- fread(file_name)
}
}
}
dat <- rbindlist(dat, idcol = "id")
dat <- separate(dat, id, c("measure", "glmname", "prewh"), sep = "__")
## calculate
## for group stats:
dat_sum <- dat %>%
group_by(prewh, measure, glmname, term, roi, subject) %>% ## average over waves
summarize(b = mean(b), .groups = "drop_last") %>%
summarize(
t_stat = t.test(b)$statistic,
p = t.test(b, alternative = "greater")$p.value,
.groups = "keep"
)
## for test-retest correlations:
dat_r <- dat %>%
pivot_wider(names_from = "wave", values_from = "b") %>% ## spread over waves
group_by(prewh, measure, glmname, term, roi) %>% ## correlate over subjects
summarize(r = cor(wave1, wave2))1 cross-run correlation
for (mod in models$crcor) {
rdm_models <-
read_model_rdm(
model = mod, measure_type = "similarity", session = session, ttype_subset = "bias"
)
p <- rdm_models %>% melt_mat %>% plot_melted_mat + labs(title = mod) + theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(p)
}
1.1 RDMs
rdm_crcor_none <- rdm$crcor__none[, , c("Vis", "SomMot", "core32", "LH_Cont_PFCl_6", "LH_Cont_Par_6"), , ]
apply(atanh(rdm_crcor_none), c("dim1", "dim2", "roi", "wave"), mean) %>%
tanh %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "group mean matrices: non-prewhitened")
rdm_crcor_none[, , , 1, ] %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "single-subject matrices: non-prewhitened")
rdm_crcor_obsall <- rdm$crcor__obsall[, , c("Vis", "SomMot", "core32", "LH_Cont_PFCl_6", "LH_Cont_Par_6"), , ]
apply(atanh(rdm_crcor_obsall), c("dim1", "dim2", "roi", "wave"), mean) %>%
tanh %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "group mean matrices: prewhitened")
rdm_crcor_obsall[, , , 1, ] %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "single-subject matrices: prewhitened")
1.2 group stats
1.2.1 non-prewhitened
dat_sum %>%
filter(measure == "crcor", prewh == "none") %>%
ggplot(aes(roi, t_stat, fill = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level t stats")
dat %>%
filter(measure == "crcor", prewh == "none") %>%
ggplot(aes(roi, b, color = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), geom = "errorbar") +
facet_grid(cols = vars(term), scales = "free") +
coord_flip() +
scale_color_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level b coefficients")
dat_sum %>%
filter(measure == "crcor", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "t_stat") %>%
ggplot(aes(lsall_1rpm, lssep_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "group-level t stats")
1.2.2 prewhitened: lsall_1rpm
dat_sum %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, t_stat, fill = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level t stats")
dat %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, b, color = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), geom = "errorbar") +
facet_grid(cols = vars(term), scales = "free") +
coord_flip() +
scale_color_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level b coefficients")
dat_sum %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "t_stat") %>%
ggplot(aes(none, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "group-level t stats")
1.3 test-retest correlations
1.3.1 non-prewhitened
dat_r %>%
filter(measure == "crcor", prewh == "none") %>%
ggplot(aes(roi, r, fill = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "crcor", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "r") %>%
ggplot(aes(lsall_1rpm, lssep_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
1.3.2 prewhitened: lsall_1rpm
dat_r %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, r, fill = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "r") %>%
ggplot(aes(none, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "crcor", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "r") %>%
ggplot(aes(obsresampbias, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
1.4 table
1.4.1 non-prewhitened
dat_sum %>%
full_join(dat_r) %>%
filter(measure == "crcor", term %in% c("conjunction", "distractor", "incongruency", "target"), prewh == "none") %>%
pivot_wider(names_from = c("glmname"), values_from = c("t_stat", "p", "r")) %>%
filter(p_lsall_1rpm < 0.01 | p_lssep_1rpm < 0.01) %>%
arrange(-t_stat_lsall_1rpm)## Joining, by = c("prewh", "measure", "glmname", "term", "roi")1.4.2 prewhitened
dat_sum %>%
full_join(dat_r) %>%
filter(measure == "crcor", term %in% c("conjunction", "distractor", "incongruency", "target"), prewh == "obsall") %>%
filter(p < 0.01) %>%
arrange(-t_stat)## Joining, by = c("prewh", "measure", "glmname", "term", "roi")All significant ROIs at \(\alpha < 0.01\) uncorrected. Sorted by t_stat_lsall_1rpm. r_* columns show test–retest correlations.
2 cross-validated euclidean
for (mod in models$cveuc) {
rdm_models <-
read_model_rdm(
model = mod, measure_type = "cvdistance", session = session, ttype_subset = "bias"
)
p <- rdm_models %>% melt_mat %>% plot_melted_mat + labs(title = mod) + theme(axis.text.x = element_text(angle = 90, hjust = 0))
print(p)
}
2.1 RDMs
rdm_cveuc_none <- rdm$cveuc__none[, , c("Vis", "SomMot", "core32", "LH_Cont_PFCl_6", "LH_Cont_Par_6"), , ]
sds <- apply(rdm_cveuc_none, c("roi", "subject", "wave"), function(x) sd(vec(x)))
rdm_cveuc_none <- sweep(rdm_cveuc_none, 3:5, sds, "/")
apply(rdm_cveuc_none, c("dim1", "dim2", "roi", "wave"), mean) %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "group mean matrices: non-prewhitened")
rdm_cveuc_none[, , , 1, ] %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "single-subject matrices: non-prewhitened")
rdm_cveuc_obsall <- rdm$cveuc__obsall[, , c("Vis", "SomMot", "core32", "LH_Cont_PFCl_6", "LH_Cont_Par_6"), , ]
sds <- apply(rdm_cveuc_obsall, c("roi", "subject", "wave"), function(x) sd(vec(x)))
rdm_cveuc_obsall <- sweep(rdm_cveuc_obsall, 3:5, sds, "/")
apply(rdm_cveuc_obsall, c("dim1", "dim2", "roi", "wave"), mean) %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "group mean matrices: prewhitened")
rdm_cveuc_obsall[, , , 1, ] %>%
melt_mat %>%
plot_melted_mat(.row = "dim1", .col = "dim2", add_greyscale_gradient = FALSE) +
scale_fill_viridis_c() +
facet_grid(vars(wave), vars(roi), scales = "free") +
theme(axis.text.x = element_text(angle = 90, hjust = 0), legend.position = "top") +
labs(title = "single-subject matrices: prewhitened")
2.2 group stats
2.2.1 non-prewhitened
dat_sum %>%
filter(measure == "cveuc", prewh == "none") %>%
ggplot(aes(roi, t_stat, fill = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level t stats")
dat %>%
filter(measure == "cveuc", prewh == "none") %>%
ggplot(aes(roi, b, color = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), geom = "errorbar") +
facet_grid(cols = vars(term), scales = "free") +
coord_flip() +
scale_color_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level b coefficients")
dat_sum %>%
filter(measure == "cveuc", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "t_stat") %>%
ggplot(aes(lsall_1rpm, lssep_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "group-level t stats")
dat_sum %>%
filter(measure == "cveuc", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "t_stat") %>%
ggplot(aes(lsall_1rpm, condition_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "group-level t stats")
2.2.2 prewhitened: lsall_1rpm
dat_sum %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, t_stat, fill = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level t stats")
dat %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, b, color = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun.data = "mean_cl_boot", position = position_dodge(width = 0.5), geom = "errorbar") +
facet_grid(cols = vars(term), scales = "free") +
coord_flip() +
scale_color_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "group-level b coefficients")
dat_sum %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "t_stat") %>%
ggplot(aes(none, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "group-level t stats")
2.3 test-retest correlations
2.3.1 non-prewhitened
dat_r %>%
filter(measure == "cveuc", prewh == "none") %>%
ggplot(aes(roi, r, fill = glmname)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "cveuc", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "r") %>%
ggplot(aes(lsall_1rpm, lssep_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "cveuc", prewh == "none") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "glmname", values_from = "r") %>%
ggplot(aes(lsall_1rpm, condition_1rpm)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
2.3.2 prewhitened: lsall_1rpm
dat_r %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
ggplot(aes(roi, r, fill = prewh)) +
geom_hline(yintercept = 0) +
stat_summary(fun = "mean", position = position_dodge(width = 0.5), geom = "col") +
facet_grid(cols = vars(term)) +
coord_flip() +
scale_fill_brewer(type = "qual", palette = 2) +
theme(legend.position = "top") +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "r") %>%
ggplot(aes(none, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
dat_r %>%
filter(measure == "cveuc", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("term", "roi"), names_from = "prewh", values_from = "r") %>%
ggplot(aes(obsresampbias, obsall)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(cols = vars(term)) +
labs(title = "test-retest correlations")
2.4 table
dat_sum %>%
full_join(dat_r) %>%
filter(measure == "cveuc", prewh == "none", term %in% c("target", "distractor", "incongruency")) %>%
pivot_wider(names_from = c("glmname"), values_from = c("t_stat", "p", "r")) %>%
filter(p_lsall_1rpm < 0.01 | p_lssep_1rpm < 0.01) %>%
arrange(-t_stat_lsall_1rpm)## Joining, by = c("prewh", "measure", "glmname", "term", "roi")dat_sum %>%
full_join(dat_r) %>%
filter(measure == "cveuc", prewh == "obsall", term %in% c("target", "distractor", "incongruency")) %>%
filter(p < 0.01) %>%
arrange(-t_stat)## Joining, by = c("prewh", "measure", "glmname", "term", "roi")All significant ROIs at \(\alpha < 0.01\) uncorrected. Sorted by t_stat_lsall_1rpm. r_* columns show test–retest correlations.
3 comparison of measures
cross-validated euclidean vs cross-run correlations
3.1 non-prewhitened euclidean
dat_sum %>%
filter(term %in% c("distractor", "incongruency", "target"), prewh == "none", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("glmname", "term", "roi"), names_from = "measure", values_from = "t_stat") %>%
ggplot(aes(crcor, cveuc)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(rows = vars(glmname), cols = vars(term)) +
labs(title = "t statistics", x = "cross-run correlation", y = "cross-validated euclidean")
dat_r %>%
filter(term %in% c("distractor", "incongruency", "target"), prewh == "none", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("glmname", "term", "roi"), names_from = "measure", values_from = "r") %>%
ggplot(aes(crcor, cveuc)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(rows = vars(glmname), cols = vars(term)) +
labs(title = "test--retest correlations", x = "cross-run correlation", y = "cross-validated euclidean")
3.2 prewhitened euclidean
dat_sum %>%
filter(term %in% c("distractor", "incongruency", "target"), prewh == "obsall", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("glmname", "term", "roi"), names_from = "measure", values_from = "t_stat") %>%
ggplot(aes(crcor, cveuc)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(rows = vars(glmname), cols = vars(term)) +
labs(title = "t statistics", x = "cross-run correlation", y = "cross-validated euclidean")
dat_r %>%
filter(term %in% c("distractor", "incongruency", "target"), prewh == "obsall", glmname == "lsall_1rpm") %>%
pivot_wider(id_cols = c("glmname", "term", "roi"), names_from = "measure", values_from = "r") %>%
ggplot(aes(crcor, cveuc)) +
geom_hline(yintercept = 0) +
geom_vline(xintercept = 0) +
geom_abline() +
geom_point(size = 0.5) +
facet_grid(rows = vars(glmname), cols = vars(term)) +
labs(title = "test--retest correlations", x = "cross-run correlation", y = "cross-validated euclidean")